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FOREWORD 


This  document  is  based  on  the  mathematical  research 
published  in  the  Technical  Report  for  Task  4 of  Contract 
No.  N6631 4-73-C-l 591  for  the  Environmental  Prediction  Research 
Facility,  Monterey,  California,  prepared  by  Philip  G.  Kesel 
and  Howard  L.  Lew it. 

The  finite  difference  equations  are  expanded  to  a 
12-level  n on - s ta gge re d model.  To  achieve  more  rapid  conver- 
gence and  to  be  more  economical  with  respect  to  computer 
storage,  the  relaxation  method  was  changed  to  sequential 
relaxation  across  vertical  planes  of  the  grid.  Stable 
solutions  were  reached  within  800  sec  on  the  Control  Data 
Corporation  (CDC)  6500  computer  in  most  cases. 
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1.  INTRODUCTION 


Initialization  of  a forecast  model  does  not  necessarily 
mean  analyzing  meteorological  data  at  given  grid  points  in 
order  to  give  a representation  of  the  actual  state  of  the 
atmosphere  as  accurately  as  possible.  It  mainly  means 
analyzing  the  data  in  such  a way  that  the  prognostic  equations 
of  the  rorecast  model  can  cope  without  producing  spurious 
gravity  waves  due  to  the  fact  that  the  initial  analysis  is  not 
in  balance  with  the  physical  constraint  impressed  on  the 
atmosphere  by  a given  set  of  prognostic  equations. 

Although  these  waves  are  eventually  suppressed  to 
reasonable  magnitudes  during  the  ti me- i terati ng  process  -- 
due  to  geostrophic  adjustment,  for  example,  in  PE  models  -- 
the  problem  remains  for  short  range  forecasts  up  to  24  hours. 

This  report  describes  a computational  concept  for 
initializing  a 12-level  PE  window  model  for  the  European- 
Mediterranean  area.  Upper-air  temperature  and  wind  analyses 
are  performed,  utilizing  the  Numerical  Variational  Analysis 
(NVA)  technique. 
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2.  APPLICATION  OF  THE  NUMERICAL 
VARIATIONAL  ANALYSIS  METHOD 


The  goal  of  the  following  upper-air  temperature  and  wind 
analyses  is  to  minimize  the  mean  square  deviation  between  the 
measured  quantities  f,  u,  and  v and  the  finally  analyzed 
quantities  T,  u,  and  v over  a 3-dimensional  domain  of  the 
atmosphere.  Thus  the  functional 

J*  =///  C ( u - u ) 2 + (v-v)2  + ( T - T ) 2 ] dxdydz  (1) 

s 

must  be  minimized. 

The  assumption  that  the  atmosphere  is  in  hydrostatic 
equilibrium  and  follows  Euler's  equation  of  motion,  and  that 
the  equation  of  state  can  be  accomplished,  leads  to  auxil- 
iary conditions  in  T,  u,  and  v which  are  introduced  into  (1) 
as  a weak  constraint  (see  Sasaki,  1970a)  giving 

j =f f f [a  ( U- U ) 2 + a ( V- V ) 2 + 3 (T-T ) 2 + yF2 

5 (2) 

+ yG2]  dxdydz 

and  its  first  variation 

6 J = ///[  a(u-u)6u  + a ( V-  V ) <5  v + 3(T-T)6T  + yF6F 

S (3) 

+ yGSG]  dxdydz. 

A priori  a,  3 and  y are  arbitrary  weights  determining  the 
degree  to  which  a term  of  the  sum  in  (2)  can  be  approximated 
by  zero  relative  to  the  others;  F and  G provide  the  dynamical 
constraint  and  are  functions  of  T,  u,  and  v. 
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2.1  DYNAMICAL  CONSTRAINT 

The  dynamical  constraints  used  in  this  case  are  derived 
from  Euler's  equation  of  motion.  Using  the  hydrostatic 
approximation,  the  equations  are  written: 

fjf  + w • wu  - fv  + = 0 (4) 


3 V 

at 


+ 


w 


\7V  + fU  + = 0 

ay 


(5) 
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p 

Introducing  the  vertical  coordinate  a = ln(— ) and  the 

P 

equation  o-f  state  into  (6),  one  gets 


i±  = R = 


9 o p 


RT. 


(7) 


Differentiating  (4)  and  (5)  with  respect  to  a and  making  use 
(7)  yields  the  thermal  wind  equations, 


9 

9 t 


/ 9 U x 9 , . 

(3^  = < 'V-'VU)  + 


3 

3a 


(fv)  - R |1 


(8) 


Jt  (£>  = - £ ( 'v-wv)  - h <fu>  - R % O) 

Assuming  that  the  left-hand  sides  of  (8)  and  (9)  are  compara- 
tively small  and  may  be  approximated  by  zero  provides  the  weak 
constraint  for  Eq.  (1).  Equations  (2)  and  (3)  are  valid  if 
we  define 
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P.  Kesel  and  H.  Lewit  (1973)  and  J.  Lewis  (1972)  used  the 
same  set  of  equations  with  the  exception  that  in  the  follow- 
ing variational  formalism,  Lewis  substituted  measured  momentum 
advection  in  (8)  and  (9). 

2.2  VARIATIONAL  TECHNIQUE  AND  ASSOCIATED  EULER  EQUATIONS 

The  necessary  and  sufficient  condition  for  minimizing 
the  functional  (2)  is  the  vanishing  of  its  first  variation 
(3). 


6J  = 0 

Substituting  (8)  and  (9)  for  F and  G in  (3)  and  assuming 
Dirichlet  conditions  at  the  boundaries,  one  obtains  after 
some  manipulation  (Kesel  and  Lewit,  1973)  the  following 
Euler  equations: 
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+ 1—  + 
3 X 
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+ r<£ 

[—  ( 
L 3X  v 

\v.\V  u)  a 

■ -2-  ( 
3y  v 

\v.\vv ) 

- £ <fv> 

9 

9 

(14) 

+ r2i^4 

ax 

4 9?T 

3y 

} = 0 

8 


a and  g are  defined  by, 


These  up  to  the  fourth  order  non-linear  set  of  partial 
differential  equations  are  solved  by  sequential  relaxation 
techniques. 


2.3  FINITE  DIFFERENCING  SCHEME  AND  RELAXATION  TECHNIQUE 

The  finite  differencing  scheme  is  set  up  for  a half-mesh 
29  X 29  grid  over  the  European-Medi terranean  area  as  shown 
in  Figure  1.  The  input  data  is  interpolated  to  half-mesh 
from  the  operational  FNWC  12-level  analysis,  the  wind  com- 
ponents being  the  geostrophic  wind  components  derived  from 
the  heigh*  fields  at  each  level.  Equations  (12)  through 
(14)  are  solved  at  each  of  the  10  interior  levels  utilizing 
the  sequential  relaxation  technique.  In  the  case  of  Eq. 

(14)  --  Helmholtz  type  of  equation  --  an  over-relaxation 
factor  of  .6  is  introduced.  Equations  (12)  and  (13)  lead  to 
the  following  general  type  of  equation  for  the  (m+l)th 
approximation  of  u and  v respectively. 


, ufj  - ♦Tt1)  + f,  fu”<)2  - (♦r:1)2)  * rt. 


i j 


i j 


1 3 


i J 


This  equation  is  linearized  giving: 


(15) 


.m+1 

i . . 

i j 


+ 2 f , 


(*Tj)2 


/ (f,  + 2 f , 


,4>m  . } 


(16) 


Equations  (12)  through  (14)  are  a simultaneous  set  of 

equations.  To  be  more  economical  with  respect  to  computer 
storage  and  computational  time,  Eqs.  (12)  through  (14)  are 
solved  on  vertical  x-o  planes  of  the  computational  grid. 
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Figure  1.  European -Mediterranean  area  (half-mesh  FNWC  grid) 


The  convergence  criteria  for  the  three  equations  are 
that  in  their  finite  difference  form,  the  respective  residuals 
are  less  than  1 at  each  grid  point.  The  simultaneous  set 
converges  toward  a stable  solution  for  T,  u,  and  v if  the 
root  mean  square  deviation  between  the  solutions  of  the  last 
3 successive  runs  through  the  atmosphere  is  less  than  .005. 

A number  of  test-runs  were  made  to  assure  that  this  condition 
is  sufficient. 


3.  SOLUTIONS 


To  ensure  that  a stable  solution  achieved  under  the 
above  mentioned  conditions  is  physically  significant,  some 
properties  of  the  solutions  have  to  be  discussed. 


3.1  CONVERGENCE  OF  THE  FINITE  DIFFERENCE  APPROXIMATION 

Since  Eqs.  (12)  through  (14)  are  a simultaneous  set  of 

non-linear  partial  differential  equations,  no  general 

convergence  criteria  could  be  found.  However,  a stable 

solution  was  achieved  with  the  finite  difference  approximation 

described  in  Section  2.3.  All  the  following  tests  have  been 

made  taking  the  FNWC  routine  analysis  for  March  16,  1972, 

0000Z , as  the  initial  guess  at  all  grid  points  interpolated 

to  half-mesh  grid  distance.  In  order  to  check  convergence 

and  the  rate  of  convergence,  the  standard  deviations  aT , a 

T \v 

for  temperature  and  wind  velocity  respectively  were  calcu- 
lated over  the  whole  domain  after  each  iteration  cycle. 


aT 


{ E 
i , j , k 


(T 


i » j , k 


(17) 


W 


= Cl 


{ E 

1 > j > k 


w 


i , j , k 


(18) 


N 25  X 25  X 10,  i and  j range  from  1 to  25,  and  k ranges 
from  1 to  10.  Figures  2 and  3 and  Table  1 show  the  results 
for  a number  of  differently  chosen  weights  a and  g . Although 
the  convergence  behavior  is  basically  similar  to  that  which 
Sasaki  (1970b)  found  investigating  shallow  water  waves,  it 
has  not  been  necessary  to  introduce  auxiliary  damping  factors 
into  the  system  in  order  to  prevent  the  iterating  process 
from  diverging.  figures  2 and  3 show  that  o-j-  and  o ^ converge 
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13 


Figure  2.  RMS  deviation  between  initial  guess  and 

NVA  temperature  field  weights  a,  6. 


W<M) 
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Figure  3.  RMS  deviation  between  initial  guess  and  NVA  wind  field  weights 


Table  1.  RMSD  between  initial  guess  and  NVA  temperature  and 
wind  fields. 
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rapidly  and  steadily  toward  a number  that  obviously  depends 
solely  on  the  weights  (a,  g)  placed  on  the  initial  data.  In 
most  cases  changes  in  Oj  and  a are  on  the  order  of  10"^ 
after  the  fifth  cycle,  except  those  cases  where  either  the 
weight  g on  the  initial  temperature  data  approaches  zero  or 
the  weight  a on  the  initial  wind  field  data  approximates  a 
value  beyond  which  convergence  can  no  longer  be  achieved. 
Another  interesting  feature  of  this  scheme  is  that  the  required 
number  of  iterations  increases  for  values  of  g approaching  0 
or  values  of  a approaching  the  above  mentioned  critical  value. 
In  addition  to  the  above,  convergence  rates  are  faster  in 
regions  where  the  smallest  differences  occur  between  the 
initial  data  and  the  final  solution. 

3.2  CONVERGENCE  TOWARD  STABLE  SOLUTIONS  WITH  RESPECT  TO 

ARBITRARILY  CHOSEN  WEIGHTS 

Sasaki  (1970c)  found  that,  in  using  the  linear  advection 
equation  as  weak  constraint  in  the  NVA , convergence  of  the 
solution  depends  on  the  ratio  between  the  weight  on  the  input 
data  and  the  weight  on  the  dynamical  constraint.  The  same 
behavior  could  be  anticipated  in  this  case,  so  several  runs 
have  been  carried  out  with  arbitrarily  chosen  weights  a and  g. 
Figure  4 shows  the  domain  in  an  a-g-plane  where  a converging 
solution  could  be  obtained;  an  "x"  marks  the  spot  where 
successful  convergence  has  been  achieved  and  no  stable  solu- 
tion could  be  found  with  weight  combinations  marked  with  0. 

It  is  interesting  to  see  that  the  value  of  g seems  to  have  no 
influence  on  con vergence, whereas  a has  to  be  chosen  larger 
than  a critical  value  to  obtain  converging  solutions.  This 
is  valid  in  the  domain  shown  in  Figure  4,  but  runs  have  been 
made  with  considerably  higher  values  of  g still  resulting 
in  convergent  solutions  with  slightly  lower  values  of  a than 
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Figure  4 . Domain  of  convergence  and  divergence 


the  critical  value  shown  in  Figure  5.  The  choices  of  the 
weights  a and  3 beyond  that  shown  in  Figure  4 left  the 
initial  fields  essentially  unchanged. 

Some  additional  characteristics  of  the  solutions  within 
the  domain  are  shown  in  Figures  5 and  6.  The  standard 
deviation  between  the  initial  guess  fields  and  the  output  NVA 
has  been  computed  according  to  Eqs.  (17)  and  (18)  with  j=3 
and  N=250 . The  isopleths  shown  were  interpolated  from  the 
data  calculated  at  the  points  x shown  in  Figure  4. 

It  is  easily  seen  from  Figure  4 that  as  the  weights  on 
the  observational  analysis  increase,  the  variation  from  this 
field  becomes  much  less  in  the  final  solution,  e.g.,  large 
values  of  the  weight  on  the  wind  field  (a)  allows  the  tempera- 
ture fields  to  vary  consi derably , whereas  the  wind  fields 
remain  essentially  unchanged. 

It  should  be  noted  that  the  RMS  deviations  shown  in 
Figures  5 and  6 are  calculated  on  a vertical  plane  next  to 
the  boundary  of  the  domain,  but  have  been  found  to  be  repre- 
sentative for  the  entire  domain.  In  order  to  show  some  of 
th°  characteristics  in  more  detail,  observational  analyses 
and  solutionsare  presented  and  discussed  in  the  next  section. 

3.3  DISCUSSION  OF  DIFFERENT  SOLUTIONS  DERIVED  FROM 

DIFFERENT  SETS  OF  WEIGHTS 

1 he  observational  analysis  should  be  described  before  the 
different  solutions  are  discussed.  The  analysis  is  from  the 
0000Z  observations  of  16  March  1972,  as  represented  by  the 
FNWC  routine  analysis  schemes  and  shown  in  Figure  7.  Figures 
8 through  11  show  the  limited  area  and  upper  air  analysis  (for 
location  of  the  systems  refer  to  Figure  1).  From  the  analysis 
it  can  be  seen  that  a European  anticyclone  stretches  to  the 
highest  tropospheric  levels;  upper-air  low-pressure  systems 
are  located  over  southwest  Europe  and  northern  Egypt;  and  the 
southeastern  portion  of  an  intensive  (945  mb)  low  pressure 
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vortex  is  present  over  the  North  Atlantic.  The  correspondi ng 
925-,  700-,  and  400-mb  temperature  and  isotach  charts  are 
labeled  "initial  guess"  in  the  sets  of  figures  shown  as 
Figures  12  through  29.  These  figures  also  show  the  NV  analysis 
for  two  different  sets  of  weights.  Without  going  into  details, 
the  following  can  be  said: 

The  technique  seems  to  straighten  out  many  of  the 
detailed  features  the  initial  guess  temperature  fields 
provide,  an  absolutely  desirable  effect  if  raw  data  assembled 
to  grid  points  is  used.  On  the  other  hand,  there  is  consider- 
able shortwave  noise  introduced  into  the  temperature  analysis 
if  the  weights  are  changed  more  and  more  in  favor  of  the 
dynamical  constraint  until  a critical  value  is  reached  and 
the  system  becomes  computationally  unstable. 

When  experimenting  with  raw  data,  it  may  turn  out  that 
a heavy  weight  on  the  dynamical  constraint  is  more  advanta- 
geous in  data-sparse  areas,  in  which  case  a smoothing 
operator  could  possibly  be  added  to  Eq.  (2)  or  run  over  the 
final  solution. 

The  isotach  charts  do  not  show  any  significant  differ- 
ences compared  to  the  initial  guess  except  at  925  mb  where 
the  velocity  gradient  belonging  to  the  North  Atlantic  low- 
pressure  system  gets  unreasonably  steep  in  the  NV  analysis 
when  the  weight  on  the  initial  data  is  decreased. 
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Figure  5.  Isopleths  of  RSMD  (°C)  between  initial  guess  and  NVA  temperature  field 
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Figure  6.  Isopleths  of  RSMD  (kt)  between  initial  guess  and  NVA  wind  field 


Figure  7.  Surface  weather  chart,  16  March  1972,  OOOOZ 
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Figure  8.  Surface  pressure,  16  March  1972,  OOOOZ 
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Figure  9.  D values  at  925  mb  (in  geopotential  decameters) 
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ligure  10.  D values  at  700  mb  (in  geopotential  decameters) 
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Figure  11.  D values  at  400  mb  (in  geopotential  decameters) 
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Figure  12.  Initial  guess  temperature  fields  (°C)  at  925  mb 
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Figure  13.  NVA  temperature  fields  (°C)  at  925  mb,  weights 
(a, 6)  = (60,25) 
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Figure  14.  NVA  temperature  fields  (°C)  at  925  mb,  weights 
(a, 3)  = (55,0) 


29 


Figure  15.  Initial  guess  temperature  fields  (°C)  at  700  mb 
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Figyre  16.  NVA  temperature  fields  (°C)  at  700mb,  weiqhts 
(0,8)  = (60,25)  * y 
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Figure  17.  NVA  temperature  fields  (°C)  at  700  mb,  weights 

(a, 3)  = (55,0) 
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Figure  18.  Initial  guess  temperature  fields  (°C)  at  400  mb 
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Figure  19.  NVA  temperature  fields  (°C)  at  400  mb,  weights, 
(a, 3)  = (60,25) 
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Figure  20.  NVA  temperature  fields  (°C)  at  400  mb,  weiqhts 
(a, 3)  = (55,0)  y 
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Figure  21.  Initial  guess  isotachs  (kt)  at  925  mb 


Figure  22.  NVA  isotachs  (kt)  at  925  mb,  weights  (a, B) =(60 ,25) 
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Figure  23.  NVA  isotachs  (kt)  at  925  mb,  weights  (a,$)=(55,0) 
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Figure  24.  Inital  guess  isotachs  (kt)  at  700  mb 
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Figure  25.  NVA  isotachs  (kt)  at  700  mb,  weights  (a, 3) =(60 , 25) 
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Figure  26.  NVA  isotachs  (kt)  at  700  mb,  weights  (a, 6) =(55,0) 
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Figure  27 . Initial  guess  isotachs  (kt)  at  400  mb 


42 


Figure  28.  NVA  isotachs  (kt)  at  400  mb,  weights  ( a, 8) = (60 , 25) 
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Figure  29.  NVA  isotachs  (kt)  at  400  mb,  weights  ( a, 3 ) = (55 ,0) 
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4.  SOME  ASPECTS  OF  THE  SPECTRAL  ANALYSIS 
OF  INITIAL  GUESS  AND  NVA  FIELDS 

Any  technique  for  analysis  of  meteorological  data  has 
to  eliminate  certain  undesirable  quantities  from  the  raw  data. 
Most  of  the  time,  with  the  techniques  in  use,  this  is 
achieved  by  applying  a low  pass  filter  to  the  raw  data  fields. 
If  the  analysis  is  to  be  used  to  initiate  a forecast  model, 
this  technique  successfully  prevents  any  major  gravity  waves 
from  building  up  during  solution  of  the  forecast  equations. 

In  use  of  the  variational  analysis  technique,  there  is 
the  possibility  that  the  technique  in  itself  provides  a low 
pass  filter  due  to  the  physical  constraint  forced  on  the  data. 

In  order  to  investigate  this,  a two-dimensional  Fourier 
analysis  was  made  from  the  NV  analysis  and  its  initial  guess. 
After  squaring  the  Fourier  coefficients,  the  ratios  between 
the  coefficients  for  the  NVA  and  the  coefficients  for  the 
initial  guess  were  calculated;  this  means  that  if  the  out- 
coming  number  was  less  than  1,  the  energy  --  with  respect 
to  that  particular  combination  of  wave  numbers  in  x and  y 
direction  --  was  reduced  in  the  NV  analysis  compared  to  the 
initial  guess.  Results  for  the  925  mb  temperature  and  wind 
fields  are  shown  in  Table  2 for  (a.  6)  = (65,25)  and  in 
Table  3 for  (a,  b)  = (55,0)  respectively. 

In  the  given  example  it  does  not  look  as  if  the  technique 
in  connection  with  the  used  physical  constraint  acts  like  a 
"filter"  for  any  specific  wavelengths.  The  only  obvious 
conclusion  that  can  be  reached  comparing  Table  2 with  Table  3 
is  that  regardless  of  the  weights  a,  b,  corrections  from  the 
initial  guess  are  applied  to  approximately  the  same  wave- 
lengths  and  the  smaller  the  weights  are,  the  more  pronounced 
are  the  corrections;  i.e.,  if  the  "enerqy"  is  reduced  with 
respect  to  a certain  wave  number,  it  will  be  reduced  even 
more  with  smaller  weights  a , b and  vice  versa. 
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Table  2.  Ratio  between  squared  Fourier  coefficients  for  NVA  over 
coefficients  for  initial  guess  for  weights  (a,  8)  = (65,25). 


* 128.6 
**  456326.2 
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Table  3.  Ratio  between  squared  Fourier  coefficients  for  the  NVA  over 
coefficients  for  initial  guess  for  weights  (a,  3)  = (55 ,0) 


* 176.9 
**  463859.1 
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It  should  be  noted  that  the  initial  guess  fields  are 
the  routine  FNWC  temperature  and  height  analyses  that  are 
most  probably  in  reasonable  agreement  with  the  physical 
constraint  used  in  this  analysis  scheme.  So  the  next  step 
was  to  add  random  noise  to  the  initial  guess.  White  noise 
with  a standard  deviation  of  1°C  and  2 knots  has  been  added 
to  the  initial  guess  temperature  and  u-  and  v-component 
fields,  respectively. 

Table  4 shows  the  ratio  between  the  squared  Fourier 
coefficients  for  the  initial  guess  fields  plus  random  noise 
over  the  squared  coefficients  for  the  initial  guess  itself. 
Table  5 shows  the  same  values,  calculated  from  the  NV  analysis 
initialized  with  fields  plus  random  noise  and  without  random 
noise. 

If  the  numbers  in  Table  5 were  equal  to  1,  it  would  mean 
that  the  NV  analysis  had  filtered  the  initially  added  random 
noise  completely.  This  is  nearly  valid  for  the  temperature 
field  and  it  can  easily  be  seen  that  a considerable  amount 
of  the  originally  added  noise  is  taken  from  u and  v component 
fields  after  application  of  the  NVA  technique. 
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Table  4.  Ratio  between  squared  Fourier  coefficients  for  initial  guess  plus 
random  noise  and  coefficients  for  initial  guess  without  random  noise. 


* 412099.8 
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Table  5.  Ratio  between  squared  Fourier  coefficients  for  NVA  initialized 
by  fields  plus  random  noise  over  coefficients  for  NVA  initialized 
by  same  fields  without  random  noise. 
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5.  CONCLUSIONS 


It  has  been  shown  that  the  NVA  technique  with  the  given 
physical  constraint  leads  to  reasonable  results.  Relatively 
small  random  noise  on  the  initial  guess  is  properly  smoothed 
However,  before  running  the  analysis  routinely,  a forecast 
model  has  to  be  initialized  in  order  to  find  the  best  set  of 
weights  on  the  initial  guess  data.  The  incorporation  of  the 
data  assembling  process  into  the  scheme  should  be  considered 
because  it  might  prove  that  otherwise  the  method  is  too  time 
consuming  if  assembling  and  analyzing  are  two  separate 
computational  processes. 
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